The averaged Kepler problem¶

We consider the Hamiltonian

$$ H(r, \theta, p_r, p_\theta) = v p_\theta + \Vert p \Vert_{g} $$

where $v$ is a constant, $p = (p_r, p_\theta)$, and $\Vert \cdot \Vert_{g}$ is the norm induced by the metric

$$ g = \mathrm{d}r^2 + m_\lambda^2(r)\, \mathrm{d}\theta^2, \quad m_\lambda^2(r) = \frac{\sin^2 r}{1 - \lambda \sin^2 r} $$

with $\lambda = 4/5$.

Along the geodesics, we have $H+p^0 = 0$. The parameter $p^0$ is constant equal to $-1$ (hyperbolic), $0$ (abnormal) or $1$ (elliptic).

Remark. We can parameterize the geodesics by the norm of the initial convector, setting $\Vert{p_0}\Vert_g = 1$. This amounts to parameterize by the initial angle $\alpha_0$: $$ p_r = \sin \alpha_0, \quad p_\theta = m_\lambda(r) \cos \alpha_0. $$ In that case, the hyperbolic geodeics are given by $$ p_\theta\, v + 1 = v\, m_\lambda(r) \cos \alpha_0 + 1 > 0. $$

Packages¶

In [1]:
import sys
sys.path.insert(1, '../')
#
import math
import numpy as np              # scientific computing tools
import wrappers                 # for compilation of Fortran Hamiltonian codes
import geometry2d.database
import geometry2d.conjugate
import geometry2d.errors
import geometry2d.geodesic
import geometry2d.plottings
import geometry2d.problem
import geometry2d.splitting
import geometry2d.utils
import geometry2d.wavefront
import nutopy as nt
import scipy
import matplotlib.pyplot as plt # for plots
%matplotlib inline
In [2]:
#
color_hyperbolic = 'red'
color_elliptic   = 'blue'
color_abnormal   = 'green'
color_cut_locus  = 'black'
color_conjugate  = 'magenta'
color_wavefront  = 'orange' 
#
color_strong_domain_2d = 'lightblue'
color_strong_domain_3d = 'black'
#
PLANE = geometry2d.plottings.Coords.PLANE
SPHERE = geometry2d.plottings.Coords.SPHERE

Fortran Hamiltonian code¶

In [3]:
# print Fortran code of the Hamiltonian: print hfun.f90 file
with open('hfun.f90', 'r') as f:
    print(f.read())

! Kepler
subroutine hfun(x, p, v, h)

    double precision, intent(in)  :: x(2), p(2), v
    double precision, intent(out) :: h

    ! local variables
    double precision :: r, th, pr, pth
    double precision :: l, m

    r   = x(1)
    th  = x(2)

    pr  = p(1)
    pth = p(2)

    l   = 4d0/5d0
    m   = sqrt(sin(r)**2/(1d0-l*sin(r)**2))

    h   = v*pth + sqrt(pr**2 + (pth/m)**2)

end subroutine hfun

Initializations¶

In [4]:
def m2(r, λ):
    return np.sin(r)**2/(1.0-λ*np.sin(r)**2)

λ_Kepler = 4.0/5.0

def g2_Kepler(r):
    return m2(r, λ_Kepler)

def initialize(v, r0):

    # Parameters
    θ0 = 0.0        # initial latitude
    t0 = 0.0        # initial time

    name = 'kepler_strong' # name of the problem

    # Initialize data
    data_file = 'data_kepler_strong.json'
    restart   = False # restart or not the computations
    data      = geometry2d.database.Data({'name': name,
                        't0': t0, 
                        'r0': r0, 
                        'θ0': θ0, 
                        'v': v}, data_file, restart)

    # Initial point
    q0 = np.array([r0, θ0])

    # Hamiltonian and derivatives up to order 3
    H = wrappers.hamiltonian(v, compile=False, display=False)

    # The Riemannian metric associated to the Zermelo problem
    # g = g1 dr^2 + g2 dθ^2
    def g(q):
        #λ  = 4.0/5.0
        r  = q[0] 
        g1 = 1.0
        g2 = g2_Kepler(r)
        return g1, g2

    # problem
    prob = geometry2d.problem.GeometryProblem2D(name, H, g, t0, q0, data, steps_for_geodesics=200)
    
    return prob
In [5]:
# function to plot the domain of strong current between two angles r1 and r2
def plot_2d_domain(fig, r1, r2):
    # from r to phi
    def phi(r):
        return r - np.pi/2.0
    ph1 = phi(r1)
    ph2 = phi(r2)
    # plot the domain: a 2d surface
    ax = fig.axes[0]
    xlims = ax.get_xlim()
    x  = np.linspace(xlims[0], xlims[1], 100)
    y1 = np.ones_like(x)*ph1
    y2 = np.ones_like(x)*ph2
    ax.fill_between(x, y1, y2, color=color_strong_domain_2d, alpha=0.5)
    return fig
In [6]:
def plot_3d_domain(fig, r1, r2, *, elevation, azimuth):
    #
    color = color_strong_domain_3d
    alpha = 0.5
    # from r to phi
    def phi(r):
        return r - np.pi/2.0
    φ1 = phi(r1)
    φ2 = phi(r2)
    # plot the domain: a 2d surface between the parallel of latitude φ1 and φ2 on the sphere
    ax = fig.axes[0]
    #
    N = 100
    u = np.linspace(0, 2 * np.pi, N)
    v = np.linspace(φ1, φ2, N)
    x = np.outer(np.cos(u), np.cos(v))
    y = np.outer(np.sin(u), np.cos(v))
    z = np.outer(np.ones(np.size(u)), np.sin(v))
    #
    x_front = np.zeros((N, N))
    y_front = np.zeros((N, N))
    z_front = np.zeros((N, N))
    x_back  = np.zeros((N, N))
    y_back  = np.zeros((N, N))
    z_back  = np.zeros((N, N))
    cam = geometry2d.plottings.get_cam(elevation, azimuth, geometry2d.plottings.dist__)
    gap = 1e-1
    for i in range(N):
        for j in range(N):
            u = np.array([x[i,j], y[i,j], z[i,j]])
            v = np.array([cam[0], cam[1], cam[2]])
            ps = np.dot(u,v)/(np.linalg.norm(u)*np.linalg.norm(v))
            #ps = x[i,j]*cam[0]+y[i,j]*cam[1]+z[i,j]*cam[2]
            
            if ps >= -gap:
                x_back[i,j] = x[i,j]
                y_back[i,j] = y[i,j]
                z_back[i,j] = z[i,j]
                x_front[i,j] = np.nan
                y_front[i,j] = np.nan
                z_front[i,j] = np.nan
            if ps <= gap:
                x_front[i,j] = x[i,j]
                y_front[i,j] = y[i,j]
                z_front[i,j] = z[i,j]
                x_back[i,j] = np.nan
                y_back[i,j] = np.nan
                z_back[i,j] = np.nan
    #
    ax.plot_surface(x_front, y_front, z_front, color=color, alpha=1) #, zorder=z_order_surface)
    ax.plot_surface(x_back, y_back, z_back, color=color, alpha=alpha) #, zorder=-z_order_surface)
    return fig

Parameterization of the geodesics¶

In [7]:
# we parameterize ||p(0)|| = 1 by the angle α0
# the norm being given by the metric g: p(0) = (sin(α0)*g1, cos(α0)*g2) = (sin(α0), cos(α0)*g2), since g1 = 1
# get the values of α0 for which h = ptheta v + 1 > 0
#
# first get the values of α for which h = ptheta v + 1 = 0
# since ptheta = m(r) cos(α0) we wil need to check if v m(r) > 1
# m(r) = sqrt(g2_Kepler(r))
def K(r, v):
    return v*np.sqrt(g2_Kepler(r))

def pθ(r, α):
    return np.cos(α)*np.sqrt(g2_Kepler(r))

def h(r, α, v):
    return pθ(r, α)*v + 1.0

def α0(r, v):
    return np.arccos(-1.0/K(r, v))

# get the values of α0 for which h > 0
def α_hyperbolique_1(r, v):
    k = K(r, v)
    if abs(k) > 1.0:
        α = α0(r, v) # between 0 and pi
        if α < np.pi/2.0:
            return α
        else:
            return 2.0*np.pi - α
    else:
        return 0.0

def α_hyperbolique_2(r, v):
    k = K(r, v)
    if abs(k) > 1.0:
        α = α0(r, v) # between 0 and pi
        if α < np.pi/2.0:
            return 2.0*np.pi - α
        else:
            return α + 2.0*np.pi
    else:
        return 2.0*np.pi

def α_hyperbolique_span(r, v, *, N=20):
    α1 = α_hyperbolique_1(r, v)
    print('α1 = ', α1)
    α2 = α_hyperbolique_2(r, v)
    print('α2 = ', α2)
    return np.linspace(α1, α2, N)

# get the values of α0 for which h < 0
def α_elliptique_1(r, v):
    k = K(r, v)
    if abs(k) > 1.0:
        α = α0(r, v) # between 0 and pi
        if α < np.pi/2.0:
            return 2.0*np.pi - α
        else:
            return α
    else:
        # raise an error: there is no elliptic domain
        raise geometry2d.errors.ArgumentValueError('There is no elliptic domain for this value of v')
    
def α_elliptique_2(r, v):
    k = K(r, v)
    if abs(k) > 1.0:
        α = α0(r, v) # between 0 and pi
        if α < np.pi/2.0:
            return α + 2.0*np.pi
        else:
            return 2.0*np.pi - α
    else:
        # raise an error: there is no elliptic domain
        raise geometry2d.errors.ArgumentValueError('There is no elliptic domain for this value of v')
    
def α_elliptique_span(r, v, *, N=20):
    α1 = α_elliptique_1(r, v)
    print('α1 = ', α1)
    α2 = α_elliptique_2(r, v)
    print('α2 = ', α2)
    return np.linspace(α1, α2, N)

def α_abnormals(r, v):
    k = K(r, v)
    if abs(k) >= 1.0:
        α = α0(r, v)
        if α < np.pi/2.0:
            return np.array([α, -α])
        else:
            return np.array([α, 2.0*np.pi-α])
    else:
        # raise an error: there is no abnormals geodesics
        raise geometry2d.errors.ArgumentValueError('There is no abnormals geodesics for this value of v')
    
def domain_strong_current(v):
    r1 = math.asin(np.sqrt(1.0/(λ_Kepler+v**2)))
    r2 = np.pi - r1
    return r1, r2

Problem definition¶

In [8]:
v    = 0.8
r0   = np.pi/2.0
prob = initialize(v, r0)
geodesic = geometry2d.geodesic.Geodesic(prob)
r1_strong, r2_strong = domain_strong_current(v) # Domain of strong current
v_Kepler = v

Geodesic flow¶

In [9]:
αspan_hyperbolique = α_hyperbolique_span(r0, v, N=21)
tf = geodesic.return_to_equator

#
fig2d = geodesic.plot(alphas=αspan_hyperbolique, tf=tf, length=1.0, view=PLANE, figsize=(5,5), color=color_hyperbolic)
ax_2d = fig2d.axes[0]
ax_2d.set_xlim(-np.pi/3, 2*np.pi + np.pi/3)
ax_2d.axvline(2*np.pi, color='k', linewidth=0.5, linestyle="dashed", zorder=geometry2d.plottings.z_order_axes)
fig2d = plot_2d_domain(fig2d, r1_strong, r2_strong)

# cameras is a list of dictionaries containing azimuth and elevetion angles
cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d = []
for cam in cameras:
    azimuth = cam['azimuth']
    elevation = cam['elevation']
    fig3d = geodesic.plot(alphas=αspan_hyperbolique, tf=tf, length=1.0, view=SPHERE, 
                          azimuth=azimuth, elevation=elevation, figsize=(3,3),
        color=color_hyperbolic)
    fig3d = plot_3d_domain(fig3d, r1_strong, r2_strong, azimuth=azimuth, elevation=elevation)
    figs3d.append(fig3d)

α1 =  4.119189204235061
α2 =  8.447181410124111

<Figure size 640x480 with 0 Axes>
In [10]:
fig2d
Out[10]:
In [11]:
figs3d[0]
Out[11]:

Abnormal geodesics¶

In [12]:
α1a, α2a = α_abnormals(r0, v)
print('α1a = ', α1a)
print('α2a = ', α2a)

α1a =  2.1639961029445254
α2a =  4.119189204235061
In [14]:
tf = geodesic.return_to_equator

#
fig2d = geodesic.plot(alphas=[α1a, α2a], tf=tf, length=1.0, view=PLANE, figsize=(5,5), color=color_abnormal)
ax_2d = fig2d.axes[0]
ax_2d.set_xlim(-np.pi/3, 2*np.pi + np.pi/3)
ax_2d.axvline(2*np.pi, color='k', linewidth=0.5, linestyle="dashed", zorder=geometry2d.plottings.z_order_axes)
fig2d = plot_2d_domain(fig2d, r1_strong, r2_strong)

# cameras is a list of dictionaries containing azimuth and elevetion angles
cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d = []
for cam in cameras:
    azimuth = cam['azimuth']
    elevation = cam['elevation']
    fig3d = geodesic.plot(alphas=[α1a, α2a], tf=tf, length=1.0, view=SPHERE, 
                          azimuth=azimuth, elevation=elevation, figsize=(3,3), color=color_abnormal)
    fig3d = plot_3d_domain(fig3d, r1_strong, r2_strong, azimuth=azimuth, elevation=elevation)
    figs3d.append(fig3d)

<Figure size 640x480 with 0 Axes>
In [15]:
fig2d
Out[15]:
In [16]:
figs3d[0]
Out[16]:

Hyperbolic geodesics¶

In [17]:
# the separatrices between hyperbolic with loops and hyperbolic without loops is given by pθ = 0

# without loops:
αs_no_loops = np.linspace(0.0, np.pi/2, 5)
αs_no_loops = αs_no_loops[[1,3]]

# with loops
αe1 = α_elliptique_1(r0, v)
αs_loops = np.linspace(np.pi/2, αe1, 7)
αs_loops = αs_loops[[1,4]]

# concatenate
αs = np.concatenate((αs_no_loops, αs_loops))

#
tf = geodesic.return_to_equator

#
fig2d = geodesic.plot(alphas=αs, tf=tf, view=PLANE, color=color_hyperbolic, figure=fig2d)

#
for fig3d in figs3d:
    geodesic.plot(alphas=αs, tf=tf, view=SPHERE, color=color_hyperbolic, figure=fig3d)
In [18]:
fig2d
Out[18]:
In [19]:
figs3d[0]
Out[19]:
In [20]:
# get the values of α0 for which the hyperbolic geodesics return to the initial point
def myfun_return_to_initial(y):
    t, α0 = y
    q = prob.geodesic(t, α0)
    return q[0] - q[-1]

t_loop_guess  = 3.0
# α0_loop_guess = middle of pi/2 and abnormal given by α1a
α0_loop_guess = 1.7 #(np.pi/2 + α1a)/2.0

y_guess = np.array([t_loop_guess, α0_loop_guess])
sol = nt.nle.solve(myfun_return_to_initial, y_guess)

t_loop_return_to_initial   = sol.x[0]
α_loop_return_to_initial_1 = sol.x[1]
α_loop_return_to_initial_2 = 2*np.pi - α_loop_return_to_initial_1

     Calls  |f(x)|                 |x|
 
         1  5.651442610204609e-02  3.448187929913334e+00
         2  3.055349141119294e-05  3.483630790700333e+00
         3  3.106629222756182e-07  3.483618968834663e+00
         4  7.495891880478653e-11  3.483618817869363e+00
         5  1.404519055154983e-14  3.483618817899530e+00

 Results of the nle solver method:

 xsol    =  [3.03970007 1.70171195]
 f(xsol) =  [-6.88338275e-15 -1.22428109e-14]
 nfev    =  5
 njev    =  1
 status  =  1
 success =  True 

 Successfully completed: relative error between two consecutive iterates is at most TolX.

Elliptic geodesics¶

In [21]:
#
αe1 = α_elliptique_1(r0, v)
αs = np.linspace(αe1, np.pi, 4)
αs = αs[1:-1]

#
tf = geodesic.return_to_equator

#
fig2d = geodesic.plot(alphas=αs, tf=tf, view=PLANE, color=color_elliptic, figure=fig2d)

#
for fig3d in figs3d:
    geodesic.plot(alphas=αs, tf=tf, view=SPHERE, color=color_elliptic, figure=fig3d)
In [22]:
fig2d
Out[22]:
In [23]:
figs3d[0]
Out[23]:

Conjugate locus¶

In [24]:
#
gap = 1e-3

# hyperbolic whitout loops
αh_no_loops_1 = -np.pi/2 + gap
αh_no_loops_2 =  np.pi/2 - gap

# hyperbolic with loops
αh_loops_1 = np.pi/2 + gap
αh_loops_2 = α_elliptique_1(r0, v) - gap
αh_loops_3 = α_elliptique_2(r0, v) + gap
αh_loops_4 = 3.0*np.pi/2 - gap

# elliptic
αe_1 = α_elliptique_1(r0, v) + gap
αe_2 = α_elliptique_2(r0, v) - gap
In [25]:
conjugate = geometry2d.conjugate.Conjugate(prob)
conj_h_no_loops   = conjugate.compute([αh_no_loops_1, αh_no_loops_2], label='hyperbolic_without_loops')
conj_h_loop_north = conjugate.compute([αh_loops_1, αh_loops_2],       label='hyperbolic_with_loops_north')
conj_h_loop_south = conjugate.compute([αh_loops_3, αh_loops_4],       label='hyperbolic_with_loops_south')
conj_elliptic     = conjugate.compute([αe_1, αe_2],                   label='elliptic')
In [26]:
# PLANE
fig2d = conjugate.plot(['hyperbolic_without_loops', 'hyperbolic_with_loops_north', 'hyperbolic_with_loops_south'],
                        view=PLANE, color=color_conjugate, figure=fig2d, linestyle='solid', linewidth=0.8)
fig2d = conjugate.plot(['elliptic'], view=PLANE, color=color_conjugate, figure=fig2d, linestyle='dashed', linewidth=0.8)

# SPHERE
for fig3d in figs3d:
    conjugate.plot(['hyperbolic_without_loops', 'hyperbolic_with_loops_north', 'hyperbolic_with_loops_south'],
                   view=SPHERE, color=color_conjugate, figure=fig3d, linestyle='solid', linewidth=0.8)
    conjugate.plot(['elliptic'], view=SPHERE, color=color_conjugate, figure=fig3d, linestyle='dashed', linewidth=0.8)
In [27]:
fig2d
Out[27]:
In [28]:
figs3d[0]
Out[28]:

Wavefronts¶

In [29]:
wavefront = geometry2d.wavefront.WaveFront(prob)
tfs = [1.0, 1.5, 2.0, 3.1, 3.5]
ws = []
for tf in tfs:
    ws.append(wavefront.compute(tf))
In [30]:
tfs = [1.0, 3.1]
for i in range(len(tfs)):
    tf = tfs[i]
    w  = ws[i]
    if i == 0:
        fig2d_wavefront = wavefront.plot(tf, view=PLANE, color=color_wavefront, figsize=(5,5))
        ax_2d = fig2d_wavefront.axes[0]
        ax_2d.set_xlim(-np.pi/2, 2*np.pi + np.pi/2)
        ax_2d.axvline(2*np.pi, color='k', linewidth=0.5, linestyle="dashed", zorder=geometry2d.plottings.z_order_axes)
        figs3d_wavefront = []
        for camera in cameras:
            azimuth = camera['azimuth']
            elevation = camera['elevation']
            fig3d_wavefront = wavefront.plot(tf, view=SPHERE, color=color_wavefront, 
                                             azimuth=azimuth, elevation=elevation, figsize=(3,3))
            figs3d_wavefront.append(fig3d_wavefront)
    else:
        wavefront.plot(tf, view=PLANE, color=color_wavefront, figure=fig2d_wavefront)
        for fig3d_wavefront in figs3d_wavefront:
            wavefront.plot(tf, view=SPHERE, color=color_wavefront, figure=fig3d_wavefront)

<Figure size 640x480 with 0 Axes>
In [31]:
fig2d_wavefront
Out[31]:
In [32]:
figs3d_wavefront[0]
Out[32]:

Cut locus¶

In [33]:
splitting = geometry2d.splitting.Splitting(prob)
In [34]:
# no loops
intersections = wavefront.self_intersections(2.0)
print(intersections['left'])
print(intersections['right'])
intersections_no_loops = intersections['right']
αspan = [0+gap, αh_no_loops_2]
print(αspan)
cut_no_loops   = splitting.compute(intersections_no_loops, αspan, label='no_loops')
t_cut_no_loops = splitting.splitting_time(cut_no_loops)

(2.0, 2.222218622020199, 4.061866589687596, array([1.57058884, 0.62521012]))
(2.0, -0.9191267057381295, 0.9203992789373718, array([1.57109573, 2.57482246]))
[0.001, 1.5697963267948967]
In [35]:
# loops: α in [α_loop_return_to_initial_1, pi/2]
intersections = wavefront.self_intersections(3.1)
print(intersections['left'])
print(intersections['right'])
intersections_loops = intersections['left']
t  = intersections_loops[0]
α1 = intersections_loops[1]
α2 = intersections_loops[2]
q  = intersections_loops[3]
intersections_loops = (t, α2, α1, q)
αspan = [α_loop_return_to_initial_1, np.pi/2+gap]
print(αspan)
cut_loops = splitting.compute(intersections_loops, αspan, label='loops')
t_cut_loops = splitting.splitting_time(cut_loops) 

(3.1, 1.6518279335268597, 4.63133677908998, array([ 1.57087215, -0.2131451 ]))
(3.1, -1.4899456366026214, 1.4900436531535355, array([1.57130931, 5.17466691]))
[1.7017119516020776, 1.5717963267948964]
In [36]:
print(α_loop_return_to_initial_2)
cut_loops.alphas

4.581473355577509
Out[36]:
array([[4.58147336, 1.70171195],
       [4.58240506, 1.70078024],
       [4.58962422, 1.69356109],
       [4.60121676, 1.68196855],
       [4.61711398, 1.66607133],
       [4.63632337, 1.64686193],
       [4.6582413 , 1.62494401],
       [4.68249347, 1.60069184],
       [4.70907165, 1.57411366],
       [4.71138898, 1.57179633]])

Optimal synthesis¶

In [37]:
# plot cut locus
fig2d_synthesis = splitting.plot(['no_loops', 'loops'], linewidth=1.0, view=PLANE, color=color_cut_locus, 
                                 figsize=(5,5))
ax_2d = fig2d_synthesis.axes[0]
ax_2d.set_xlim(-np.pi/2, 2*np.pi+np.pi/2)
ax_2d.axvline(2*np.pi, color='k', linewidth=0.5, linestyle="dashed", zorder=geometry2d.plottings.z_order_axes)
fig2d_synthesis = plot_2d_domain(fig2d_synthesis, r1_strong, r2_strong)

cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d_synthesis = []
for cam in cameras:
    azimuth = cam['azimuth']
    elevation = cam['elevation']
    fig3d_synthesis = splitting.plot(['no_loops', 'loops'], linewidth=1.0, view=SPHERE, color=color_cut_locus, 
                                     azimuth=azimuth, elevation=elevation, figsize=(3,3))
    fig3d_synthesis = plot_3d_domain(fig3d_synthesis, r1_strong, r2_strong, azimuth=azimuth, elevation=elevation)
    figs3d_synthesis.append(fig3d_synthesis)

<Figure size 640x480 with 0 Axes>
In [38]:
# add hyperbolic conjugate points
fig2d_synthesis = conjugate.plot(['hyperbolic_with_loops_north', 'hyperbolic_with_loops_south', 
                                  'hyperbolic_without_loops'], linewidth=1.0, view=PLANE, 
                                 color=color_conjugate, figure=fig2d_synthesis, linestyle='solid')

for fig3d_synthesis in figs3d_synthesis:
    conjugate.plot(['hyperbolic_with_loops_north', 'hyperbolic_with_loops_south', 
                    'hyperbolic_without_loops'], linewidth=1.0, view=SPHERE, 
                   color=color_conjugate, figure=fig3d_synthesis, linestyle='solid')
In [39]:
# add abnormal geodesics until the cusp point defined by the boundary of the domain of strong current
α1a, α2a = α_abnormals(r0, v)
#
print('α1a = ', α1a)
print('α2a = ', α2a)
print('r1_strong = ', r1_strong)
print('r2_strong = ', r2_strong)

# the abnormal associated to α1a meets the boundary r = r2_strong
def myfun(t):
    return prob.geodesic(t, α1a)[-1][0] - (r2_strong - 1e-5)
opt = nt.nle.Options(Display='off')
t_cusp_α1a = nt.nle.solve(myfun, 1.0, options=opt).x

# the abnormal associated to α2a meets the boundary r = r1_strong
def myfun(t):
    return prob.geodesic(t, α2a)[-1][0] - (r1_strong + 1e-5)
t_cusp_α2a = nt.nle.solve(myfun, 1.0, options=opt).x

# t_cusp_abnormal
def t_cusp_abnormal(α):
    if α == α1a:
        return t_cusp_α1a
    elif α == α2a:
        return t_cusp_α2a
    else:
        raise geometry2d.errors.ArgumentValueError('α must be equal to α1a or α2a')

# plot
fig2d_synthesis = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=1.0,
                                color=color_abnormal, figure=fig2d_synthesis)
fig2d_synthesis = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=1.0,
                                color=color_cut_locus, linestyle='dashed', figure=fig2d_synthesis)

for fig3d_synthesis in figs3d_synthesis:
    geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=1.0,
                  color=color_abnormal, figure=fig3d_synthesis)
    geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=1.0,
                    color=color_cut_locus, linestyle='dashed', figure=fig3d_synthesis)

α1a =  2.1639961029445254
α2a =  4.119189204235061
r1_strong =  0.9851107833377455
r2_strong =  2.1564818702520476
In [40]:
fig2d_synthesis
Out[40]:
In [41]:
figs3d_synthesis[0]
Out[41]:
In [42]:
restart_path_loop_1 = False
if prob.data.contains('path_loop_1') and not restart_path_loop_1:
    path_loop_1 = prob.data.get('path_loop_1')
    yout_loop_1 = np.array(path_loop_1['yout'])
    αout_loop_1 = np.array(path_loop_1['αout'])
    print('path_loop_1 loaded from data')
else:
    # add hyperbolic geodesics with loops until the cut point defined by the intersection with the abnormals
    # for α in [pi/2, αh_loops_1] the intersection is with the abnormal α1a
    def myfun_itersection_1(y, α):
        th, ta, r, θ = y
        q  = np.array([r, θ])
        eq = np.zeros(4)
        qf = prob.geodesic(th, α[0])[-1]
        eq[0:2] = q - qf
        eq[2:4] = q - prob.geodesic(ta, α1a)[-1]
        return eq

    # we start with α close to α1a since we know a good approximation of the intersection
    α0_ = α1a - 1e-3
    th = t_cusp_α1a
    ta = t_cusp_α1a
    q  = prob.geodesic(th, α0_)[-1]
    r  = q[0]
    θ  = q[1]

    y_guess = np.array([th, ta, r, θ])
    sol = nt.nle.solve(lambda y : myfun_itersection_1(y, [α0_]), y_guess)
    y0 = sol.x

    print('th = ', y0[0])
    print('ta = ', y0[1])
    print('t_cusp_α1a = ', t_cusp_α1a)
    sol = nt.path.solve(myfun_itersection_1, y0, [α0_], [α_loop_return_to_initial_1])
    yout_loop_1, αout_loop_1 = sol.xout, sol.parsout
    path_loop_1 = {'yout': yout_loop_1.tolist(), 'αout': αout_loop_1.tolist()}
    prob.data.update({'path_loop_1': path_loop_1})
    print('path_loop_1 saved in data')

path_loop_1 loaded from data
In [43]:
restart_path_loop_2 = False
if prob.data.contains('path_loop_2') and not restart_path_loop_2:
    path_loop_2 = prob.data.get('path_loop_2')
    yout_loop_2 = np.array(path_loop_2['yout'])
    αout_loop_2 = np.array(path_loop_2['αout'])
    print('path_loop_2 loaded from data')
else:
    # add hyperbolic geodesics with loops until the cut point defined by the intersection with the abnormals
    # for α in [pi/2, αh_loops_1] the intersection is with the abnormal α1a
    def myfun_itersection_2(y, α):
        th, ta, r, θ = y
        q  = np.array([r, θ])
        eq = np.zeros(4)
        qf = prob.geodesic(th, α[0])[-1]
        eq[0:2] = q - qf
        eq[2:4] = q - prob.geodesic(ta, α2a)[-1]
        return eq

    # we start with α close to α1a since we know a good approximation of the intersection
    α0_ = α2a + 1e-3
    th = t_cusp_α2a
    ta = t_cusp_α2a
    q  = prob.geodesic(th, α0_)[-1]
    r  = q[0]
    θ  = q[1]

    y_guess = np.array([th, ta, r, θ])
    sol = nt.nle.solve(lambda y : myfun_itersection_2(y, [α0_]), y_guess)
    y0 = sol.x

    print('th = ', y0[0])
    print('ta = ', y0[1])
    print('t_cusp_α2a = ', t_cusp_α2a)
    sol = nt.path.solve(myfun_itersection_2, y0, [α0_], [α_loop_return_to_initial_2])
    yout_loop_2, αout_loop_2 = sol.xout, sol.parsout
    path_loop_2 = {'yout': yout_loop_2.tolist(), 'αout': αout_loop_2.tolist()}
    prob.data.update({'path_loop_2': path_loop_2})
    print('path_loop_2 saved in data')

path_loop_2 loaded from data
In [44]:
# add hyperbolic geodesics with loops until the cut point defined by the intersection with the abnormals
def t_meet_abnormal(α):
    α = α % (2.0*np.pi)
    if α <= α1a and α >= α_loop_return_to_initial_1:
        times  = yout_loop_1[:, 0]
        alphas = αout_loop_1[:, 0] % (2*np.pi)
    elif α >= α2a and α <= α_loop_return_to_initial_2:
        times  = yout_loop_2[:, 0]
        alphas = αout_loop_2[:, 0]  % (2*np.pi)
    else:
        raise geometry2d.errors.ArgumentValueError('α must be in [α_loop_return_to_initial_1, α1a] U [α2a, α_loop_return_to_initial_2]')
    t_meet = scipy.interpolate.interp1d(alphas, times, kind='cubic')(α)
    t_equator = geodesic.first_return_to_equator(α)
    return min(t_meet, t_equator)

# plot
αspan1 = np.linspace(α_loop_return_to_initial_1+1e-3, α1a-1e-3, 4)
αspan1 = αspan1[0:-1]
αspan2 = np.linspace(α_loop_return_to_initial_2-1e-3, α2a+1e-3, 4)
αspan2 = αspan2[0:-1]
αspan  = np.concatenate((αspan1, αspan2))

fig2d_synthesis = geodesic.plot(alphas=αspan, tf=t_meet_abnormal, view=PLANE, linewidth=0.5,
                                color=color_hyperbolic, figure=fig2d_synthesis)

for fig3d_synthesis in figs3d_synthesis:
    geodesic.plot(alphas=αspan, tf=t_meet_abnormal, view=SPHERE, linewidth=0.5,
                  color=color_hyperbolic, figure=fig3d_synthesis)
In [45]:
fig2d_synthesis
Out[45]:
In [46]:
# add hyperbolic geodesics without loops until the cut locus
αspan = np.linspace(αh_no_loops_1, αh_no_loops_2, 11)

fig2d_synthesis = geodesic.plot(alphas=αspan, tf=t_cut_no_loops, view=PLANE, linewidth=0.5,
                                color=color_hyperbolic, figure=fig2d_synthesis)

for fig3d_synthesis in figs3d_synthesis:
    geodesic.plot(alphas=αspan, tf=t_cut_no_loops, view=SPHERE, linewidth=0.5,
                  color=color_hyperbolic, figure=fig3d_synthesis)
In [47]:
fig2d_synthesis
Out[47]:
In [48]:
# add hyperbolic geodesics with loops but not intersection with abnormals until the cut locus
αspan1 = np.linspace(np.pi/2+1e-3, α_loop_return_to_initial_1, 3)
αspan2 = np.linspace(α_loop_return_to_initial_2+1e-3, 3.0*np.pi/2-1e-3, 3)
αspan  = np.concatenate((αspan1, αspan2))

fig2d_synthesis = geodesic.plot(alphas=αspan, tf=t_cut_loops, view=PLANE, linewidth=0.5,
                                color=color_hyperbolic, figure=fig2d_synthesis)

for fig3d_synthesis in figs3d_synthesis:
    geodesic.plot(alphas=αspan, tf=t_cut_loops, view=SPHERE, linewidth=0.5,
                  color=color_hyperbolic, figure=fig3d_synthesis)
In [49]:
fig2d_synthesis
Out[49]:
In [50]:
figs3d_synthesis[0]
Out[50]:

Level lines of the Hamiltonian¶

In [51]:
gap = 1e-3
delta = 0.025
rs  = np.arange(gap, np.pi-gap, delta)
prs = np.arange(-1.0, 1.0, 1*delta)
RS, PRS = np.meshgrid(rs, prs)
PTS = np.zeros((len(prs), len(rs)))

def pthetafun(r, pr, pzero, i, j):
    t = 0.0
    q = np.array([r, 0.0]) # theta is a cyclic variable
    def fun(ptheta):
        p = np.array([pr, ptheta])    
        return prob.Hamiltonian(t, q, p) + pzero
    # solve fun(ptheta) = 0
    # for ptheta_guess, use a value close to PTS[i, j]: do not access PTS[i, j] directly or out of bounds
    m, n = PTS.shape
    if i == 0:
        i_ = 1
    elif i == m - 1:
        i_ = m - 2
    else:
        i_ = i+1
    if j == 0:
        j_ = 1
    elif j == n - 1:
        j_ = n - 20
    else:
        j_ = j+1
    pheta_guess = PTS[i_, j_]
    #ptheta_sol = scipy.optimize.fsolve(fun, pheta_guess)
    options = nt.nle.Options(Display='off')
    sol = nt.nle.solve(fun, pheta_guess, options=options)
    return sol.x

pzero = -1.0
for j, r in enumerate(rs):
    for i, pr in enumerate(prs):
        H = pthetafun(r, pr, pzero, i, j)
        PTS[i,j] = H
    
fig, ax = plt.subplots()
CS = ax.contour(RS, PRS, PTS)
#ax.clabel(CS, inline=True, fontsize=10)
#ax.set_title('Simplest default with labels')
pass
In [52]:
delta = 0.01
prmax = 3.0
gap = 1e-3
rs  = np.arange(0+gap, np.pi-gap, delta)
prs = np.arange(-prmax, prmax, delta/prmax)
RS, PRS = np.meshgrid(rs, prs)

def pthetafun(r, pr, pzero, sign):
    m2_ = g2_Kepler(r)
    v_  = v_Kepler
    a = (1/m2_) - v_**2
    b = -2*v*pzero
    c = pr**2 - pzero**2
    delta = b**2 - 4*a*c
    if delta < 0.0:
        return np.nan
    else:
        ptheta1 = (-b + np.sqrt(delta))/(2*a)
        ptheta2 = (-b - np.sqrt(delta))/(2*a)
        if -pzero - v_ * ptheta1 >= 0.0 and -pzero - v_ * ptheta2 >= 0.0:
            if sign == 1:
                m = max(ptheta1, ptheta2)
                if m >= 0:
                    return m
                else:
                    return np.nan
            elif sign == -1:
                m = min(ptheta1, ptheta2)
                if m <= 0:
                    return m
                else:
                    return np.nan
            else:
                raise geometry2d.errors.ArgumentValueError('sign must be equal to 1 or -1')
        elif -pzero - v_ * ptheta1 >= 0.0:
            m = ptheta1
            if sign * m >= 0:
                return m
            else:
                return np.nan
        elif -pzero - v_ * ptheta2 >= 0.0:
            m = ptheta2
            if sign * m >= 0:
                return m
            else:
                return np.nan
        else:
            return np.nan
 
# hyperbolic
pzero = -1.0
#
PTS_HYP_NEG = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
    for i, pr in enumerate(prs):
        pθ = pthetafun(r, pr, pzero, -1)
        PTS_HYP_NEG[i,j] = pθ

PTS_HYP_POS = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
    for i, pr in enumerate(prs):
        pθ = pthetafun(r, pr, pzero, 1)
        PTS_HYP_POS[i,j] = pθ

# elliptic
pzero = 1.0
PTS_ELL = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
    for i, pr in enumerate(prs):
        pθ = pthetafun(r, pr, pzero, -1)
        PTS_ELL[i,j] = pθ
    
# abnormal
pzero = 0.0
PTS_ABN = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
    for i, pr in enumerate(prs):
        pθ = pthetafun(r, pr, pzero, -1)
        PTS_ABN[i,j] = pθ
In [53]:
restart_calcul_level_lines = True
if prob.data.contains('level_lines') and not restart_calcul_level_lines:
    level_lines = prob.data.get('level_lines')
    PTS_HYP_NEG = np.array(level_lines['PTS_HYP_NEG'])
    PTS_HYP_POS = np.array(level_lines['PTS_HYP_POS'])
    PTS_ELL = np.array(level_lines['PTS_ELL'])
    PTS_ABN = np.array(level_lines['PTS_ABN'])
    print('level_lines loaded from data')
else:
    level_lines = {'PTS_HYP_NEG': PTS_HYP_NEG.tolist(), 'PTS_HYP_POS': PTS_HYP_POS.tolist(),
                     'PTS_ELL': PTS_ELL.tolist(), 'PTS_ABN': PTS_ABN.tolist()}
    prob.data.update({'level_lines': level_lines})
    print('level_lines saved in data')

level_lines saved in data
In [54]:
# function to plot the domain of strong current between two angles r1 and r2
def plot_2d_domain_level_lines(ax, r1, r2):
    # plot the domain: a 2d surface
    ylims = ax.get_ylim()
    x  = np.linspace(r1, r2, 100) 
    y1 = np.ones_like(x)*ylims[0]
    y2 = np.ones_like(x)*ylims[1]
    ax.fill_between(x, y1, y2, color=color_strong_domain_2d, alpha=0.5)
    return ax
In [58]:
#
fig, ax = plt.subplots()

#
#CS = ax.contour(RS, PRS, PTS_HYP_NEG, [-4, -3, -2, -1, -0.25], colors=color_hyperbolic, linewidths=1)
#ax.clabel(CS, inline=True, fontsize=10)
CS = ax.contour(RS, PRS, PTS_HYP_POS, [0.02, 0.25, 0.5, 0.75], colors=color_hyperbolic, linewidths=1, linestyles='dashed')
#ax.clabel(CS, inline=True, fontsize=10)

# add geodesics
r0 = np.pi/2.0
θ0 = 0.0
tf = 6.12
pzero = -1.0
pr0s = np.linspace(1.25, 3.0, 4)
pθ0s = []
for pr0 in pr0s:

    pθ0 = pthetafun(r0, pr0, pzero, -1)
    pθ0s += [pθ0]
    q0 = np.array([r0, θ0])
    p0 = np.array([pr0, pθ0])
    t0 = 0.0
    N = 200
    tspan = list(np.linspace(t0, tf, N+1))
    qs, ps = prob.extremal(t0, q0, p0, tspan)
    r  = [q[0] for q in qs]
    pr = [p[0] for p in ps]
    ax.plot(r, pr, color=color_hyperbolic, linewidth=1.0, linestyle='solid')

#
ax.set_xlim(0.0, np.pi)
ax.set_ylim(-3.0, 3.0)
plot_2d_domain_level_lines(ax, r1_strong, r2_strong)

#
print('pθ0s = ', pθ0s)

#
pass

pθ0s =  [-0.3228913274337241, -1.1266351099358802, -1.9642005576625727, -2.8172904669025316]
In [93]:
#
fig, ax = plt.subplots()

#
CS = ax.contour(RS, PRS, PTS_ELL, [-20, -8, -6, -4, -3], colors=color_elliptic, linewidths=1, linestyles='solid')
plot_2d_domain_level_lines(ax, r1_strong, r2_strong)

if False:
    
    # add geodesics
    θ0 = 0.0
    pr0 = 0.0
    tf = 5
    pzero = 1.0
    r0s = np.linspace(1.1, 1.5, 4)
    pθ0s = []

    for r0 in r0s:

        pθ0 = pthetafun(r0, pr0, pzero, -1)
        pθ0s += [pθ0]
        q0 = np.array([r0, θ0])
        p0 = np.array([pr0, pθ0])
        print('q0 = ', q0)
        print('p0 = ', p0)
        t0 = 0.0
        N = 200
        tspan = list(np.linspace(t0, tf, N+1))
        qs, ps = prob.extremal(t0, q0, p0, tspan)
        r  = [q[0] for q in qs]
        pr = [p[0] for p in ps]
        ax.plot(r, pr, color=color_hyperbolic, linewidth=1.0, linestyle='solid')

#
pass
In [96]:
#
fig, ax = plt.subplots()

#
CS = ax.contour(RS, PRS, PTS_ABN, [-6, -4, -3, -2, -1, -0.25], colors=color_abnormal, linewidths=1, linestyles='solid')
plot_2d_domain_level_lines(ax, r1_strong, r2_strong)

if False:
    
    # add geodesics
    θ0 = 0.0
    r0 = np.pi/2.0
    tf = 5
    pzero = 0.0
    pr0s = np.linspace(0.1, 2.5, 4)
    pθ0s = []

    for pr0 in pr0s:

        pθ0 = pthetafun(r0, pr0, pzero, -1)
        pθ0s += [pθ0]
        q0 = np.array([r0, θ0])
        p0 = np.array([pr0, pθ0])
        print('q0 = ', q0)
        print('p0 = ', p0)
        t0 = 0.0
        N = 200
        tspan = list(np.linspace(t0, tf, N+1))
        qs, ps = prob.extremal(t0, q0, p0, tspan)
        r  = [q[0] for q in qs]
        pr = [p[0] for p in ps]
        ax.plot(r, pr, color=color_hyperbolic, linewidth=1.0, linestyle='solid')

#
pass
In [84]:
# 3D plots
from matplotlib import cm
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
surf = ax.plot_surface(RS, PRS, PTS_HYP_POS, cmap=cm.coolwarm, linewidth=0, antialiased=False)
#ax.set_zlim(-20.0, 0.0)
In [71]:
min(PTS_HYP_NEG[1000,:])
Out[71]:
-202.83241871800467